####################################################################
# AGER FIRE GENERATOR
####################################################################
# Created: 11.1.2011
# Last Update: 5.17.2012 (Update: Add Oak Knoll Fuel Moisture Info)
####################################################################
    ##METHODS:
          # Read MC1 monthly ERC values
          # Generate daily ERC values from monthly values using std dev using B. Johnson methods
          # Determine the probability of a fire from Preisler logistic regression
          # Determine the size of the fire from Greenfell splines and exp distribution
          # Determine wind speed and wind direction from Wind Rose Gust data
          # Convert the fire size to burn period based on FConstMTT calibrations
          # Determine fuel moisture values from regression, and generate fuel moisture files
          # Create fire list file
          # Create fire list file with probability switch

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#  INFORMATION YOU NEED TO RUN PROGRAM:
#   1) Set filename path (working directory); place all files in this location (line 61)
#   2) Install necessary packages (lines 57-58)
#   3) Files: a) monthly MC1 input data file (vcfutercgpred.csv)
#             b) historical fire events file (erc-raws-ignitions-binary2.csv)
#             c) historical fire ignitions file (erc-avgfiresize.csv)
#             d) hitorical wind gust file (BrushCreek_WindRose.csv)
#             e) calibration burn period and windspeed data (calfirelist250.csv)
#             f) historical fuel moisture file (villagecrk_fuelmoisture.csv)
#   4) Start and end dates for future MC1 monthly data
#   5) Set location to output fuel moisture files, see line 659
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#   CHOOSE YOUR CLIMATE SCENARIO (**search on keyword 'climate'**)
#---------------------------------------------------------------------
# You must specify the climate scenario in the following locations:
#   1) Line 78: select model by specifying column number
#   2) Line 119: select model by specifying column number
#   3) Line 659: specify path for fuel moisture files for Randig
#   4) Line 703: change the firelist filename based on model
#   5) Line 727: change the firelist ignition filename based on model
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#############################################
#  CHOOSE CRAN MIRROR
#--------------------------------------------
# Either select the CRAN mirror manually via your GUI interface or,
# run code below!

chooseCRANmirror(graphics=getOption("menu.graphics"))   ## RUN ONLY ONCE ##

######################################
#  PACKAGES NEEDED:
######################################
#   1. DATA.TABLE
######################################

install.packages("data.table")   ## RUN ONLY ONCE ##
install.packages("lattice")      ## RUN ONLY ONCE ##

#Set working directory
setwd("/Users/cevers/Dropbox/cnh/fire/FireGenerator/")

#Load libraries
library (data.table)
library (lattice)

##################################################
# GENERATE DAILY MC1 ERC VALUES FROM MONTHLY MC1 ##############################
##################################################
  #To generate daily values within R....
  #Generate list of all years, months and days for our data record
  #For each day calculate the variance in ERC

mc1.month.mat <- read.csv("vcfutercgpred.csv") #Read monthly ERC values

## ---> Specify climate model!

erc.pred.mc1 <- mc1.month.mat[,4] #Define ERC variable where col4 = HADLEY, col5 = CSIRO, col6 = MIROC
yr.mc1 <- mc1.month.mat[,3] #Define year
mo.mc1 <- mc1.month.mat[,2] #Define month

#------------------------------------------
#Create data set with all years, months and days......

#Enter start and end dates for sequence
dateseq <- seq(as.Date("2007/1/1"), as.Date("2099/12/31"), by="day")

jday <- strftime(dateseq, format="%j") #Create julian day-of-year variable
month <- strftime(dateseq, format="%m") #Define month variable
month <- as.integer(month)
day <- strftime(dateseq, format="%d") #Define day variable
day <- as.integer(day)
year <- strftime(dateseq, format="%Y") #Define year variable (with century)
year <- as.integer(year)
fireday <- strftime(dateseq, format="%B-%d") #Define fireday variable

final.dates <- cbind(dateseq, year, month, day, jday)
final.dates <- as.data.frame(final.dates)
#---------------------------------------------

#Merge final dates file with original monthly ERC file
mc1.month.mat2 <- mc1.month.mat[,-1] #First delete id variable
mc1.daily.mat <- merge(final.dates, mc1.month.mat2, by=c("year", "month"), all=TRUE)
mc1.obs <- nrow(mc1.daily.mat) #Get number of observations
#There are 33,968 observations in the dataset

#R converts year and month to factor during the merge step- these must be converted back to maintain the correct
#vector order
mc1.daily.mat$year.int <- as.numeric(levels(mc1.daily.mat$year)[mc1.daily.mat$year]) #Convert factor to integer!
mc1.daily.mat$mon.int <- as.numeric(levels(mc1.daily.mat$month)[mc1.daily.mat$month]) #Convert factor to integer!
mc1.daily.mat.final <- mc1.daily.mat[,-1] #Remove old year factor variable
mc1.daily.mat.final <- mc1.daily.mat.final[,-1] #Remove old month factor variable

#Resort matrix to original order by year, month and day
mc1.daily.sort <- mc1.daily.mat.final[order(mc1.daily.mat.final$year.int, mc1.daily.mat.final$mon.int, mc1.daily.mat.final$day) , ]

## ---> Specify climate model!

erc.pred.mc1 <- mc1.daily.sort[,4] #Define ERC variable again where col4 = HADLEY, col5 = CSIRO, col6 = MIROC

#Calculate standard deviation of ERC (Equation from B. Johnson)
erc.sd <- -1.80293+3.22763*log(erc.pred.mc1+0.1)

#Calculate daily ERC (Method from B. Johnson)
#set.seed(7) #Inactivated to develop variability and multiple list files
erc.daily <- erc.pred.mc1 + (erc.sd*rnorm(mc1.obs)) #Specify # of observations

#Replace negative values with zero and round
erc.mc1 <- round((ifelse (erc.daily<0, 0, erc.daily)), digits=4)

#Create dataframe with year, month, day, julian day of year, and predicted daily ERC
mc1.daily.input <- cbind(year, month, day, jday, erc.mc1)
mc1.daily.input <- as.data.frame(mc1.daily.input)

##NOTE: Final ERC vector = erc.mc1
##NOTE: Final input dataframe = mc1.daily.input

############################
# DETERIMINE FIRE IGNITION ###################################
############################
  #Based on ERC for 4 RAWS stations and FPA Ignition data for ignition study area
  #RAWS Stations (North to South):
        #Village Creek
        #Mt. Yoncalla
        #Silver Butte
        #Provolt Seed Orchard
  #There are 2,184 missing dates between 8/5/1988 and 12/31/2008, thus only 5,270 days in record
  #Data set contains ERC values for a given RAWS Station Zone matched spatially to ignitions
  #Thus dates repeat for each station (i.e., 21,080 days )
  #Data set contains date, RAWS Station (or region), historical daily ERC value, regional mean ERC and fire event

erc.mat <- read.csv("erc-raws-ignitions-binary2.csv") #Historical ignition data

erc.hist <- erc.mat[,3] #Define historical ERC variable
region <- erc.mat[,5] #Define RAWS station zones
event <- erc.mat[,6] #Define event 0 = no fire, 1 = fire
mean.ERC <- erc.mat[,7] #Define regional mean ERC

erc.hist.std <- erc.hist-mean.ERC #Standardize ERCs based on regional mean

  #Logistic Regression Model with Region - fire events and historical ERCs
    
erc.logit <- glm(event ~ as.numeric(erc.hist.std)+region-1, family = binomial("logit"))
summary(erc.logit)

#-------------------------------------------------------------------
#  CREATE PROBABILITY
#-------------------------------------------------------------------
#  log[p0/(1-p0)] -b x ERCbar  + b x ERC
#
#     p0 = study area ignitions/all ignitions
#     b = slope from logistic regression (coefficient of erc.hist)
#       = 0.052185 (for these data)
#     ERCbar = historical mean ERC from study area
#     ERC = erc.mc1
#
# Equations created by Haiganoush Preisler 2.16.12
#-------------------------------------------------------------------

beta <- coefficients(erc.logit) #Create matrix of coefficients from glm output
b <- subset(beta, beta > 0, select = beta) #Select erc.hist coefficient

new.logit <- (log(0.013/(1-0.013))-(b*27)+(b*erc.mc1)) #Where b = 0.052185
P <- exp(new.logit)/(1+exp(new.logit)) #Calculate probability

#------------------------------------------------------------------------
#   STUDY AREA INFORMATION
#------------------------------------------------------------------------
#   Ignition Area: 
#           3,879,338 hectares 
#           4,222 ignitions, 201 ignitions/year
#           Average # of days per/year with >=1 ignition = 66 (35-160)
#   Lebanon Area: 
#           101,887 hectares 
#           2.6% of ignition area
#           45 ignitions, 2.1 ignitions/year
#           Fraction of observed fires in entire ignitions area = 0.012
#           Mean historical ERC = 34 (for ignitions)
#   Eugene Area: 
#           82,100 hectares 
#           2.1% of ignition area 
#           56 ignitions, 2.6 ignitions/year
#           Fraction of observed fires in entire ignition area = 0.013
#           Mean historical ERC = 32 (for ignitions)
#   Village Creek Historical Mean ERC = 27 (excluding zeros)
#   Mt. Yoncalla Historical Mean ERC = 24 (excluding zeros)
#   Silver Butte Historical Mean ERC = 27 (excluding zeros)
#   Provolt Seed Orchard Historical Mean ERC = 29 (excluding zeros)
#------------------------------------------------------------------------

pred.fire.mc1 <- round(P, digits=4) #Redefine probability variable

hist(erc.mc1) #Histogram of MC1 daily ERCs
hist(erc.hist) #Histogram of historical daily ERCs

plot (erc.mc1, pred.fire.mc1)

fire.prob <- cbind(erc.mc1, pred.fire.mc1) #Create matrix of mc1 data and fire probs

  ##NOTE: Final fire probability vector = pred.fire.mc1
  ##NOTE: Final fire probability matrix = fire.prob

########################
# DETERIMINE FIRE SIZE ############################
########################
  #Code written by Isaac Greenfell
  #Establish relationship between ERC and Fire Size from Historical Data

hist.mat <- read.csv("erc-avgfiresize.csv") #Historical FPA ignition data

ercseq <- hist.mat[,3] #Define ERC variable
fsizeseq <- hist.mat[,7] #Define Fire Size variable

#Convert fire size in acres to hectares
fsizeseq <- fsizeseq * 0.404685642

par(mfrow=c(1,1))

## SMOOTH SPLINE - NonParametric Regression Model
  #Changing df will adjust how wiggly the line is
  #Relationship between Fire Size and ERC
sm1 <- smooth.spline(fsizeseq~ercseq, df = 5)

# Domain over which fit is to be made
predseq <- 1:120 

# Calculate fitted values over range of ERC
pred.sizeseq <- predict(sm1, predseq)

plot(ercseq, fsizeseq)
lines(pred.sizeseq$x, pred.sizeseq$y, col="green")

#Assuming there is a fire using MC1 derived daily ERCs... 
mean.simseq <- predict(sm1, erc.mc1) #Calculate fitted values over range of MC1 ERC values
plot(mean.simseq$x, mean.simseq$y)

#Generate random number from exponential distribution with rate
  #This is conditional on day
out.sim <- rexp(erc.mc1, 1/mean.simseq$y)

#Plot new simulation data versus historical
par(mfrow = c(2,1))

max.ha <- 800

plot(mean.simseq$x, out.sim, ylim = c(0, max.ha))
plot(ercseq, fsizeseq, ylim = c(0, max.ha))

fire.size <- round(out.sim, digits=2) #Redefine variable

  ##NOTE: Final fire size vector (ha) = fire.size

##########################################################
# GENERATE WIND SPEED AND AZIMUTH FROM WIND ROSE GUST   #####################
##########################################################
#Based on Brush Creek RAWS Station historical data
#Years 2000 to 2010, units = mph
#Data are daily summary time series

wind.rose <- read.csv("BrushCreek_WindRose.csv") #Read wind data
rose.obs <- nrow(wind.rose)

w.date <- as.Date((wind.rose[,1]), "%m/%d/%Y") #Define date
w.month <- strftime(w.date, format="%m") #Define month variable
w.year <- as.integer(strftime(w.date, format="%Y")) #Define year variable (with century)
w.myr <- strftime(w.date, format="%m-%Y") #Create Month-by-Year variable

wind.rose$w.myr <- w.myr #Add month-by-year variable to dataframe

w.rose <- as.data.frame(wind.rose[, c(7,8)]) #Select needed variables

w.rose.mean <- aggregate(w.rose, list(as.factor(w.rose$w.myr)), mean, na.rm = T) #Calculate mean wind gust by month-year
w.rose.mean$w.mo <- as.integer(substr(w.rose.mean$Group.1, 1, 2)) #Create month variable

w.rose.sm <- as.data.frame(w.rose.mean[, c(2, 4)]) #Create dataframe with only needed variables; NULLs CAUSE WARNINGS
w.rose.sm[10,1] <- 24 #Replace missing values for one year with average of other years -January
w.rose.sm[21,1] <- 26.68 #Replace missing velus for one year with average of other years -February

#Create separate arrays for each month
w.jan <- subset(w.rose.sm, w.mo == 1)
  w.jan <- as.vector(w.jan[,1])
w.feb <- subset(w.rose.sm, w.mo == 2)
  w.feb <- as.vector(w.feb[,1])
w.mar <- subset(w.rose.sm, w.mo == 3)
  w.mar <- as.vector(w.mar[,1])
w.apr <- subset(w.rose.sm, w.mo == 4)
  w.apr <- as.vector(w.apr[,1])
w.may <- subset(w.rose.sm, w.mo == 5)
  w.may <- as.vector(w.may[,1])
w.jun <- subset(w.rose.sm, w.mo == 6)
  w.jun <- as.vector(w.jun[,1])
w.jul <- subset(w.rose.sm, w.mo == 7)
  w.jul <- as.vector(w.jul[,1])
w.aug <- subset(w.rose.sm, w.mo == 8)
  w.aug <- as.vector(w.aug[,1])
w.sep <- subset(w.rose.sm, w.mo == 9)
  w.sep <- as.vector(w.sep[,1])
w.oct <- subset(w.rose.sm, w.mo == 10)
  w.oct <- as.vector(w.oct[,1])
w.nov <- subset(w.rose.sm, w.mo == 11)
  w.nov <- as.vector(w.nov[,1])
w.dec <- subset(w.rose.sm, w.mo == 12)
  w.dec <- as.vector(w.dec[,1])

gust <- rep(1, mc1.obs) #Define variable

#Loop through each observation to randomly select max. wind gust by month
for (i in 1:mc1.obs){
  
  if (month[i] == 1) {
    gust[i] <- sample(w.jan, size=1);
  } else if (month[i] == 2) {
    gust[i] <- sample(w.feb, size=1);
  } else if (month[i] == 3) {
    gust[i] <- sample(w.mar, size=1);
  } else if (month[i] == 4) {
    gust[i] <- sample(w.apr, size=1);
  } else if (month[i] == 5) {
    gust[i] <- sample(w.may, size=1);
  } else if (month[i] == 6) {
    gust[i] <- sample(w.jun, size=1);
  } else if (month[i] == 7) {
    gust[i] <- sample(w.jul, size=1);
  } else if (month[i] == 8) {
    gust[i] <- sample(w.aug, size=1);
  } else if (month[i] == 9) {
    gust[i] <- sample(w.sep, size=1);
  } else if (month[i] == 10) {
    gust[i] <- sample(w.oct, size=1);
  } else if (month[i] == 11) {
    gust[i] <- sample(w.nov, size=1);
  } else if (month[i] == 12) {
    gust[i] <- sample(w.dec, size=1)
  }
} # ends loop
                
w.spd <- round(gust, digits=0) #Redefine wind gust as w.spd

#Wind Direction

dir.quad <- rep(0, rose.obs)
wind.rose$wind.dir[is.na(wind.rose$wind.dir)] <- 999
                                                      
wind.dir <- wind.rose$wind.dir
for (i in 1:rose.obs){
  
  if (wind.dir[i] <= 60) {
    dir.quad[i] <- 1;
  } else if ((wind.dir[i] > 60) & (wind.dir[i] <= 120)) {
    dir.quad[i] <- 2;
  } else if ((wind.dir[i] > 120) & (wind.dir[i] <= 180)) {
    dir.quad[i] <- 3;
  } else if ((wind.dir[i] > 180) & (wind.dir[i] <= 240)) {
    dir.quad[i] <- 4;
  } else if ((wind.dir[i] > 240) & (wind.dir[i] <= 300)) {
    dir.quad[i] <- 5
  } else if ((wind.dir[i] > 300) & (wind.dir[i] <= 360)) {
    dir.quad[i] <- 6;
  } else if (wind.dir[i] == 999) {
    dir.quad[i] <- 9 
  }
}

#Determine most frequent quadrat for each month
freq.table <- table(w.month, dir.quad)
freq.table

#RESULTS (month followed by most freq quadrant):
#THIS IS HARDWIRED IN- IF DATA ARE CHANGED THEN THE UPCOMING LOOP WILL NEED TO BE FIXED!
#jan 4, feb 4, march 4, apr 4, may 6, jun 6, jul 6, aug 6, sep 6, oct 4, nov 4, dec 4

azimuth <- rep(NA, mc1.obs) #Define azimuth variable

#Create loop that generates random number of quadrant compass range based on above results
for (i in 1:mc1.obs){
  
  #For Months BTW May and September, Quadrant is 6 (compass 301-360)
  if ((month[i] > 4) & (month[i] < 10)) {
    azimuth[i] <- sample((seq(301:360)), size=1);
  #For All Other Months, Quadrant is 4 (compass 181-204)
  } else {
    azimuth[i] <- sample((seq(181:240)), size=1);
    }
} #ends loop

  ##NOTE: Final wind speed vector = w.spd
  ##NOTE: Final wind direction vector = azimuth


##############################
# RUN BURN PERIOD TRANSLATOR #####################
##############################
  # Calibration data created from Randig-FConstMTT simulations, 250 ignition multiplier
    # landscape file used: eidu4_90mTS0000.lcp created on 3.19.12
    # resolution: 90; crown fire method: ScottRheinhart

#Data set with burn period, windspeed, and fire size (ha)

calibrate <- read.csv("calfirelist250.csv") #Read randig calibration data

fm80 <- calibrate[grep("80", calibrate$FMFile),] #Create subsetted dataframe for fuelmoisture 80
fm30 <- calibrate[grep("30", calibrate$FMFile),] #Create subsetted dataframe for fuelmoisture 30

#Create Windspeed-by-BurnPeriod variable to average data:
fm80$by <- paste(as.character(fm80$WindSpeed), as.character(fm80$BurnPer), sep="_") 
fm80.sub <- fm80[, c(4,5,10,11)] #Select only needed variables

#Calculate average firesize by Windspeed-BurnPeriod
fm80.avg <- aggregate(fm80.sub, list(as.factor(fm80.sub$by)), mean)

#------------------------------------------
#Graph Fire Size by Burn Period to examine relationships
#NOTE: WINDSPEED DOES NOT IMPACT BURN PERIOD IN THIS STUDY AREA

attach(fm80.avg) #Grab variable names
par(mfrow = c(1,1)) #Set graph output parameters

xyplot(Hectares ~ BurnPer, data=fm80.avg, groups=WindSpeed, cex = 1.5, pch = c(15,17,6,7,8,15,17,16), 
       col=c("red", "darkgrey", "green", "orange", "deepskyblue2", "darkred", "purple", "blue4"), key = 
         list(title = "WindSpeed", text=list(c("5","10", "15", "20", "25", "30", "35","40")),
              x = .05, y=.95, corner = c(0,1), 
              border = FALSE, lines = FALSE, 
              points=list(pch=c(15,17,6,7,8,15,17,16), 
                          col=c("red", "darkgrey", "green", "orange", "deepskyblue2", "darkred", "purple", "blue4")))
       )
#------------------------------------------

#Generate new means based on burn period only!
fm80.bp <- aggregate(fm80.sub, list(as.factor(fm80.sub$BurnPer)), mean)

attach(fm80.bp) #Grab variable names
par(mfrow = c(1,1)) #Set graph output parameters
plot(BurnPer ~ Hectares)

bp.reg <- lm((BurnPer~Hectares+I(Hectares^2))) #Polynomial linear regression model
summary(bp.reg)

sizeha <- fire.size #Redefine fire.size data

preddat <- data.frame(Hectares = sizeha) #Create new dataframe with redefined fire size

pred.bp <- predict(bp.reg, newdata = preddat) #Predict new burn period with calibration data

plot(preddat[,1], pred.bp) #Plot fire size by predicted burn period

burnp <- round(pred.bp, digits=0) #Redefine burnperiod variable and round

  ##NOTE: Final burn period vector = burnp

####################################
# GENERATE FUEL MOISTURE FILE LINK ##########################
####################################
#Village Creek ERC-fuel moisture data from Charles McHugh
    # villagecreek_frisk_file_May_Oct.txt
#Oak Knoll ERC-fuel moisture data from Charles McHugh
    # oak_knoll_frisk_file_may_oct.txt
#IMPORTANT NOTE: In order to capture high historical ERCs, Oak Knoll data will be used to predict
  # ERC Values above 60!!

fuel <- read.csv('villagecrk_fuelmoisture.csv')#Village Creek Fuel Moisture

erc.fuel <- fuel[,1] #Define ERC variable
fm1 <- fuel[,2] #Define 1 hr fuel moisture variable
fm10 <- fuel[,3] #Define 10 hr fuel moisture variable
fm100 <- fuel[,4] #Define 100 hr fuel moisture variable
fm1000 <- fuel[,5] #Define 1000 hr fuel moisture variable
fmherb <- fuel[,6] #Define herbaceous fuel moisture variable
fmwood <- fuel[,7] #Define woody fuel moisture variable

fuel.oak <- read.csv('oakknoll_fuelmoisture.csv')#Oak Knoll Fuel Moisture
fuel.high <- fuel.oak[fuel.oak$ERC>=60, ]#Grab only data where ERC >=60

erc.fuel.oak <- fuel.high[,1] #Define ERC variable -- Oak Knoll
fm1.oak <- fuel.high[,2] #Define 1 hr fuel moisture variable -- Oak Knoll
fm10.oak <- fuel.high[,3] #Define 10 hr fuel moisture variable -- Oak Knoll
fm100.oak <- fuel.high[,4] #Define 100 hr fuel moisture variable -- Oak Knoll
fm1000.oak <- fuel.high[,5] #Define 1000 hr fuel moisture variable -- Oak Knoll
fmherb.oak <- fuel.high[,6] #Define herbaceous fuel moisture variable -- Oak Knoll
fmwood.oak <- fuel.high[,7] #Define woody fuel moisture variable -- Oak Knoll

## FUEL MOISTURE - 1 HR FUELS------------------------
fm1.reg <- lm((fm1~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fm1.reg)

fm1.reg.oak <- lm((fm1.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model-- Oak Knoll
summary(fm1.reg.oak)

par(mfrow = c(4,1)) #Set graph output parameters
max.erc <- 120 #Set x-axis maximum value

plot(erc.fuel, fm1, xlim = c(0, max.erc), main = "Village Creek 1hr Fuel Moisture") #Plot Village Creek
lines(erc.fuel, predict(fm1.reg), col="blue") 

plot(erc.fuel.oak, fm1.oak, xlim = c(60, max.erc), main = "Oak Knoll Fuel Moisture- High ERC") #Plot Oak Knoll
lines(erc.fuel.oak, predict(fm1.reg.oak), col="blue")

tempx <- erc.mc1 #Redefine erc.mc1 data
hightempx <- subset(erc.mc1, erc.mc1 > 59.9999) #Redefine erc.mc1 data -- Oak Knoll

newdat <- data.frame(erc.fuel = tempx) #Create new dataframe with redefined mc1
newdat.high <- data.frame(erc.fuel.oak = hightempx) #Create new dataframe with redefined mc1 -- Oak Knoll

pred.fm1 <- predict(fm1.reg, newdata = newdat) #Predict new 1hr fuel moistures with erc.mc1 data
pred.fm1.oak <- predict(fm1.reg.oak, newdata=newdat.high) #Predict new 1hr fuel moistures with erc.mc1 data --  Oak Knoll

plot(newdat[,1], pred.fm1, main = "Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 1hr fuel moisture
plot(newdat.high[,1], xlim = c(60, max.erc), pred.fm1.oak, main = "Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 1hr fuel moisture

## FUEL MOISTURE - 10 HR FUELS------------------------------------------------------------------------
fm10.reg <- lm((fm10~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fm10.reg)

fm10.reg.oak <- lm((fm10.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model
summary(fm10.reg.oak)

plot(erc.fuel, fm10, xlim = c(0, max.erc), main = "Village Creek 10hr Fuel Moisture")#Plot Village Creek
lines(erc.fuel, predict(fm10.reg), col="blue")

plot(erc.fuel.oak, fm10.oak, xlim = c(60, max.erc), main = "Oak Knoll Fuel Moisture- High ERC")#Plot Village Creek
lines(erc.fuel.oak, predict(fm10.reg.oak), col="blue")

pred.fm10 <- predict(fm10.reg, newdata = newdat) #Predict new 10hr fuel moistures with erc.mc1 data
pred.fm10.oak <- predict(fm10.reg.oak, newdata=newdat.high)#Predict new 10hr fuel moisture with erc.mc1 data -- Oak Knoll

plot(newdat[,1], pred.fm10, main = "Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 10hr fuel moisture
plot(newdat.high[,1], xlim = c(60, max.erc), pred.fm10.oak, main = "Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 10hr fuel moisture - Oak Knoll

#Fuel moistures above 2.7 are incorrectly predicted!  These data will be replaced.
#Fuel moistures should reach asymptote and not continue to drop.
#Replace values where fuel moisture <2.7 with 2.69 based on historical Oak Knoll data

pred.fm10.oak2 <- ifelse(pred.fm10.oak <2.7, 2.69, pred.fm10.oak) #Correct asymptote

plot(newdat.high[,1], xlim = c(60, max.erc), pred.fm10.oak2, main = "Predicted Fuel Moisture for High ERC - corrected",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 10hr fuel moisture - Oak Knoll

## FUEL MOISTURE - 100 HR FUELS------------------------------------------------------------------------
fm100.reg <- lm((fm100~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fm100.reg)

fm100.reg.oak <- lm((fm100.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model-- Oak Knoll
summary(fm100.reg.oak)

plot(erc.fuel, fm100, xlim = c(0, max.erc), main="Village Creek 100hr Fuel Moisture")#Plot Vilage Creek
lines(erc.fuel, predict(fm100.reg), col="blue")

plot(erc.fuel.oak, fm100.oak, xlim = c(60, max.erc), main="Oak Knoll Fuel Moisture- High ERC")#Plot Vilage Creek
lines(erc.fuel.oak, predict(fm100.reg.oak), col="blue")

pred.fm100 <- predict(fm100.reg, newdata = newdat) #Predict new 100hr fuel moistures with erc.mc1 data
pred.fm100.oak <- predict(fm100.reg.oak, newdata = newdat.high) #Predict new 100hr fuel moistures with erc.mc1 data-- Oak Knoll

plot(newdat[,1], pred.fm100, main="Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 100hr fuel moisture
plot(newdat.high[,1], pred.fm100.oak, xlim = c(60, max.erc), main="Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 100hr fuel moisture -- Oak Knoll

## FUEL MOISTURE - 1000 HR FUELS------------------------------------------------------------------------
#Note that 1000 hr fuels are NOT included in the fuel moisture files

## FUEL MOISTURE - HERBACEOUS FUELS------------------------------------------------------------------------
fmherb.reg <- lm((fmherb~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fmherb.reg)

fmherb.reg.oak <- lm((fmherb.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model
summary(fmherb.reg.oak)

plot(erc.fuel, fmherb, xlim = c(0, max.erc), main="Village Creek Herb. Fuel Moisture")
lines(erc.fuel, predict(fmherb.reg), col="blue")

plot(erc.fuel.oak, fmherb.oak, xlim = c(60, max.erc), main="Oak Knoll Fuel Moisture- High ERC")
lines(erc.fuel.oak, predict(fmherb.reg.oak), col="blue")

pred.fmherb <- predict(fmherb.reg, newdata = newdat) #Predict new herb. fuel moistures with erc.mc1 data
pred.fmherb.oak <- predict(fmherb.reg.oak, newdata=newdat.high)

plot(newdat[,1], pred.fmherb, main="Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted herb. fuel moisture
plot(newdat.high[,1], xlim = c(60, max.erc), pred.fmherb.oak, main="Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted")

#Fuel moistures above 80 are incorrectly predicted!  These data will be replaced.
#Fuel moistures should reach asymptote.
#Replace fuel moisture values where ERC >=80 with 2 based on historical Oak Knoll data
pred.fmherb.oak2 <- ifelse(hightempx>=80, 2, pred.fmherb.oak)
plot(newdat.high[,1], pred.fmherb.oak2) #Plot mc1 ERC by corrected fuel moisture

## FUEL MOISTURE - WOODY FUELS------------------------------------------------------------------------
fmwood.reg <- lm((fmwood~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fmwood.reg)

fmwood.reg.oak <- lm((fmwood.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model
summary(fmwood.reg.oak)

plot(erc.fuel, fmwood, xlim = c(0, max.erc), main="Village Creek Woody Fuel Moisture")
lines(erc.fuel, predict(fmwood.reg), col="blue")

plot(erc.fuel.oak, fmwood.oak, xlim = c(60, max.erc), main="Oak Knoll Fuel Moisture- High ERC")
lines(erc.fuel.oak, predict(fmwood.reg.oak), col="blue")

pred.fmwood <- predict(fmwood.reg, newdata = newdat) #Predict new woody fuel moistures with erc.mc1 data
plot(newdat[,1], pred.fmwood, xlim = c(60, max.erc), main="Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted woody fuel moisture

pred.fmwood.oak <- predict(fmwood.reg.oak, newdata = newdat.high) #Predict new woody fuel moistures with erc.mc1 data
plot(newdat.high[,1], xlim = c(60, max.erc), pred.fmwood.oak, main="Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted") #Plot mc1 ERC by predicted woody fuel moisture

#Fuel moistures above 87 are incorrectly predicted!  These data will be replaced.
#Fuel moistures should reach asymptote.
#Replace fuel moisture values where ERC >=87 with 60 based on historical Oak Knoll data
pred.fmwood.oak2 <- ifelse(hightempx>=87, 60, pred.fmwood.oak)
plot(newdat.high[,1], pred.fmwood.oak2) #Plot mc1 ERC by corrected fuel moisture
                        
#------------------------------------------------------------------

#Create dataframe with erc.mc1 and all fuel categories
newfuel <- cbind(erc.mc1, pred.fm1, pred.fm10, pred.fm100, pred.fmherb, pred.fmwood)
newfuel<- as.data.frame(newfuel)
lowfuel <- subset(newfuel, newfuel$erc.mc1 < 60)

highfuel <- as.data.frame(cbind(erc.mc1=hightempx, pred.fm1=pred.fm1.oak, pred.fm10=pred.fm10.oak2,
                                pred.fm100=pred.fm100.oak, pred.fmherb=pred.fmherb.oak2, 
                                pred.fmwood = pred.fmwood.oak2))

##Calculate mean value for each fuel moisture category by ERC value, and then round
fuel.mean1 <- round(aggregate(lowfuel, list(round(lowfuel$erc.mc1, digits=0)), mean), digits=0)
fuel.mean2 <- round(aggregate(highfuel, list(round(highfuel$erc.mc1, digits=0)), mean), digits=0)

##ACTIVATE THE FOLLOWING TWO LINES OF CODE TO VIEW FUEL MOISTURE VALUES BY ERC!!!
fuel.mean <- rbind(fuel.mean1, fuel.mean2) #Create table of fuel moisture values
#write.table (fuel.mean, file="fuelmoisture_hadley.txt", sep=",", quote= F, row.names=FALSE) 

#Set the location to output the files:
## ---> Specify climate model!

path <- "C:\\runfiregen\\randig\\HAD\\"

filect <- nrow(fuel.mean)

#Create fuel moisture files, one file per ERC value
for (i in 1:filect){
  
  ID <- fuel.mean$erc.mc1[i]
  
  data.now <- paste(("0"), " ",
              paste(fuel.mean$pred.fm1[i]), " ",
              paste(fuel.mean$pred.fm10[i]), " ",
              paste(fuel.mean$pred.fm100[i]), " ",
              paste(fuel.mean$pred.fmherb[i]), " ",
              paste(fuel.mean$pred.fmwood[i]))
  
  write(data.now, paste(path, paste("FlamMapFuelMoistures", ID,".fms", sep=""), sep=""))
  
} # ends loop

#Assign fuel moisture file to each daily ERC value (data row)
fuel.file <- rep(NA, mc1.obs)

erc.mc1.int <- round(erc.mc1, digits = 0)#Create ERC as integer
fmlabel <- as.data.frame(erc.mc1.int)#Create label for fuel moisture file

for (i in 1: mc1.obs){
  
  fuel.file[i] <- paste("FlamMapFuelMoistures", fmlabel$erc.mc1.int[i], ".fms", sep="")  
} #end loop

###########################
# GENERATE FIRE LIST FILE ##########################
###########################
  #Create table that outputs year, julian day of year,
    #  fire probability, burn period, wind speed, azimuth and fuel moisture file 
    #  for each day in MC1 input file

## year, prob, julianday ,burn period, wind (mph), azimuth, FMS File, daily ERC
fire.list = data.table (yr=year-2007, prob = pred.fire.mc1, julian=jday, burnper = burnp, windspeed = w.spd,
                       azimuth = azimuth, FMfile = fuel.file, ERC = erc.mc1, size = fire.size)

## ---> Specify climate model!

write.table (fire.list, file="firelist_HAD1.txt", sep=",", quote= F, row.names=FALSE) #For Envision

################################################
# GENERATE FIRE LIST FILE WITH IGNITION SWITCH ##########################
################################################
#Create table that outputs year, julian day of year,
#  fire probability, burn period, wind speed, azimuth and fuel moisture file 
#  for each day in MC1 input file after implementing ignition swtich

#Determine ignition! ##Ignition is determined in Envision- 
#RUN ONLY WHEN YOU NEED PROBABILITY

rand.prob <- (runif(mc1.obs, 0, 1)) #Generate random number
ignition <- ifelse ((pred.fire.mc1 >= rand.prob), 1, 0)

#Create data frame of fileid (chardate), erc, ignition, windspeed, azimuth, burn period and fuel file
fire.list.ig <- as.data.frame(cbind(yr=year-2007, prob=pred.fire.mc1, julian=jday, burnper=burnp, windspeed=w.spd, 
                                      azimuth, FMfile=fuel.file, ERC=erc.mc1,size= fire.size, ignition))

#Grab only ignitions
fires.only <- as.data.frame(fire.list.ig[fire.list.ig$ignition==1, ])

## ---> Specify climate model!

write.table (fires.only, file="firelist_HAD1_ig.txt", sep=",", quote= F, row.names=FALSE) #For Envision
